function bpdf = custom_betapdf(x, a, b)
% custom_betapdf  - PDF of the beta distribution
%
% FORMAT:       bcdf = custom_betacdf(x, a, b)
%
% Input fields:
%
%       x           values for which the beta PDF is computed
%       a, b        beta function parameters
%
% Output fields:
%
%       bpdf        beta pDF of for (x, a, b)

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 3 || ...
   ~isnumeric(x) || ...
   ~isnumeric(a) || ...
   ~isnumeric(b) || ...
    isempty(a) || ...
    isempty(b)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(x)
    bpdf = [];
    return;
end
osize = size(x);
if numel(x) == 1
    if numel(a) > 1
        osize = size(a);
    elseif numel(b) > 1
        osize = size(b);
    end
end
if ~isa(x, 'double')
    x = double(x(:));
else
    x = x(:);
end
if ~isa(a, 'double')
    a = double(a(:));
else
    a = a(:);
end
if ~isa(b, 'double')
    b = double(b(:));
else
    b = b(:);
end

% build output
bpdf = zeros(osize);
bpdf(~(a > 0) | ~(b > 0) | isnan(x)) = NaN;

% regular cases
k = find(x > 0 & x < 1 & a > 0 & b > 0);
if ~isempty(k)
    if numel(a) == 1 && ...
        numel(b) == 1
        bpdf(k) = exp((a - 1) .* log(x(k)) + ...
		    (b - 1) .* log(1 - x(k))) ./ beta(a, b);
    elseif numel(a) == 1
        bpdf(k) = exp((a - 1) .* log(x(k)) + ...
		    (b(k) - 1) .* log(1 - x(k))) ./ beta(a, b(k));
    elseif numel(b) == 1
        bpdf(k) = exp((a(k) - 1) .* log(x(k)) + ...
		    (b - 1) .* log(1 - x(k))) ./ beta(a(k), b);
    else
        bpdf(k) = exp((a(k) - 1) .* log(x(k)) + ...
		    (b(k) - 1) .* log(1 - x(k))) ./ beta(a(k), b(k));
    end
end

% reshape
bpdf = reshape(bpdf, osize);
